重要な注意¶
- このNotebookを開いただけの状態では，編集結果は保存されないので，各自，「ファイル」メニューから「ドライブにコピーを保存」してください．
- 操作方法の詳細は Python / Colab 超入門 で勉強してください．
１．使用するライブラリの読込¶
In [1]:
#この枠をクリックしてアクティブにしてから，Shiftを押しながらEnterを押すと，枠内のコードが実行されます．以下同じです．
# '#'記号のあとの内容を「コメントアウト」と呼び，実行時に無視されます．
import numpy as np #数値計算ライブラリ
from math import sin, cos #低機能だが計算が速い三角関数
from scipy.integrate import odeint #常微分方程式ライブラリ
import matplotlib.pyplot as plt #グラフ作成ライブラリ
from matplotlib.animation import FuncAnimation #アニメーションライブラリ
from matplotlib import rc #グラフ調整ライブラリ
from IPython.display import HTML #ユーザーインターフェース拡張ライブラリ
２．運動方程式の数値解を求める関数¶
パラメータの設定¶
In [2]:
M = 2/3 #台車の質量
m = 1/3 #振り子の質量
l = 1 #振り子の長さ
数値解を求める関数¶
In [3]:
def Solve(xx0, force):
### 運動方程式を1階化した微分方程式の定義
def eom(xx,t):
#注 python は配列の添字が0から
x = xx[0] #台車の変位
dx = xx[1] #台車の速度
th = xx[2] #振り子の角度
dth = xx[3] #振り子の角速度
ft = force(xx) #制御力
A = np.array( #行列
[
[M+m, m*l*cos(th)],
[np.cos(th), l]
]
)
bb = np.array( #右辺のベクトル
[m*l*(dth**2)*sin(th)+ft, 9.8*sin(th)]
)
A_inv = np.linalg.inv(A) #逆行列
hh = np.dot(A_inv, bb) #逆行列とベクトルの積
dxx = np.array( #1階微分のベクトル
[dx, hh[0], dth, hh[1]]
)
return dxx
###ここまで
### 時間軸を表す等差数列
n = 200 #時刻の数
ts = np.linspace(0, 25, n) #0秒から25秒までn等分
### 差分解法（数値積分）
xxs = odeint(eom, xx0, ts)
return (ts,xxs)
３．実際に数値解を求める¶
外力を与える関数（ひとまず外力は０にしておく）¶
In [4]:
def force(xx):
return 0
初期値の設定（ちょっと傾いて静止）¶
In [5]:
xx0 = np.array([0, 0, 0.5, 0]) #時刻t=0 のx, dx, th, dth
そのときの数値解を求める¶
In [6]:
ts, xxs = Solve(xx0, force)
ベクトル xx := (x, dx, th, dth) を **状態ベクトル** という．
In [7]:
print(xxs) #各行が状態ベクトル，下に向かって時間が進行
４．数値解の振動波形を確認する¶
振動波形を描く関数¶
In [8]:
def Plot_Wave(ts, xxs):
plt.plot(ts,xxs[:,0],label='x') #台車の変位
plt.plot(ts,xxs[:,1],label='dx/dt') #台車の速度
plt.plot(ts,xxs[:,2],label='th') #振り子の角度
plt.plot(ts,xxs[:,3],label='d(th)/dt') #振り子の角速度
plt.xlabel('t')
plt.ylabel('States')
plt.legend()
数値解の振動波形¶
In [9]:
Plot_Wave(ts, xxs)
５．数値解からアニメーションを作る¶
数値解をアニメーションする関数¶
In [10]:
def Animate(xxs, title='xxxxxxX'):
### アニメーション用のグラフ用紙を用意する
fig, ax = plt.subplots(figsize=(8,3)) #グラフ用紙を作る
plt.close() #ひとまず表示OFF
ax.set_xlim(-4,4) #グラフの縦軸の範囲
ax.set_ylim(-1.5,1.5) #グラフの横軸の範囲
ax.grid() #グリッドon
ax.set_xlabel('X') #横軸のラベル
ax.set_ylabel('Y') #縦軸のラベル
ax.set_title(title) #タイトル（学籍番号）
line1, = ax.plot([], [], lw=2) #空の描画
### i 行目の状態ベクトルで描画データを更新する関数
def update_anim(i):
x = xxs[i,0] #さっき計算したxxsの1列目（台車変位，添字は0）
th = xxs[i,2] #さっき計算したxxsの3列目（振子角度，添字は2）
#振り子支点の位置ベクトル
XM = np.array([x, 0])
#振り子先端の位置ベクトル
Xm = XM + l*np.array([sin(th), cos(th)])
#振り子の線分を描画
line1.set_data([XM[0],Xm[0]],[XM[1],Xm[1]]) #座標データの更新
### アニメーションデータの作成
n = len(xxs) #データの行数（時間きざみの総数）
anim = FuncAnimation(fig, update_anim, interval=100, frames=n)
### アニメーションの出力形式の設定
rc('animation', html='jshtml')
return anim
実際にアニメーションする¶
- 表示されるまで少し計算時間がかかる．
- 表示されたら，再生ボタンを押すと，アニメーションが始まる．
In [11]:
Animate(xxs)
Out[11]:
- こうしたシミュレーションでは，運動方程式に書いてない効果は，当然，再現されない．
- 重力に抗する力は与えていないので，振り子は落下し，それに連動して台車も動く．
- 床を作成していないので，振り子は下まで落ちて，スイングする．
- 摩擦や空気抵抗がないので，スイングは持続する．
実習12.1（改）¶
- 最初から順に実行して，ここまでくれば，実習12.1は完了しています．
実習12.2（改）¶
- 次のコードセルの
forceとxx0を，
def force(xx):
return 20*xx[2]
xx0 = np.array([0, -3, 0.5, 1.2])のように書き換えて実行せよ．
In [12]:
def force(xx):
return 0
xx0 = np.array([0, 0, 0.5, 0])
ts, xxs = Solve(xx0, force)
Animate(xxs)
Out[12]:
実習12.3（改）¶
- 次のコードセルの
forceとxx0を，
def force(xx):
return 20*xx[2] + 2*xx[3]
xx0 = np.array([0, -3, 0.3, 1.2])のように書き換えて実行せよ．
In [13]:
def force(xx):
return 0
xx0 = np.array([0, 0, 0.5, 0])
ts, xxs = Solve(xx0, force)
Animate(xxs)
Out[13]:
実習12.4（改）¶
次のコードセルの
KとLの数値を自分で書き換え（必要ならxx0も），次の2種類の運動を実現せよ．(1) 振動しながら立位に向う運動．
(2) 振動しないで立位に向う運動．
In [14]:
def force(xx):
K = 20
L = 2
return K*xx[2] + L*xx[3]
xx0 = np.array([0, -3, 0.3, 1.2])
ts, xxs = Solve(xx0, force)
Animate(xxs)
Out[14]: